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Abstract 

Starting from the balance equations of mass, momentum and energy we formulate 
an integral ID model for a poly-disperse mixture injected in the atmosphere. We write 
all the equations, either in their most general formulation or in the more simplihed, 
taking particular care in considering all the underlying hypothesis in order to make clear 
when it is possible and appropriate to use them. Moreover, we put all the equations 
in a non-dimensional form, making explicit all the dimensionless parameters that drive 
the dynamics of these phenomena. In particular, we hnd parameters to measure: the 
goodness of the Boussinesq approximation, the injected mass flow, the column stability 
and his eventual collapse, and the importance of the atmospheric stratification, the initial 
kinetic energy and the gravitational potential energy. We show that setting to zero some 
of these parameters, it is possible to recover some of the existing jet and plume models 
for single-phase flows. Moreover, we write a simplihed set of equations for which it is 
possible to hnd analytical solutions that can be used to describe also the dynamics of 
multiphase “weak-plumes”. 

Starting from the paper HH the study on jets and plumes has been carried out by 
a lot of different researcher involved in a variety of disciplines. Indeed, these kind of 
phenomena are quite ubiquitous in nature... 


1 The main assumptions. 

In order to use the Dusty Gas model we have to assume: 

• Local equilibrium. 

• All the phases, either solid or gaseous, move with the same velocity held 
u{x,t). [S] shows that this assumption it is valid even for the solid phase if 
the Stokes time Tg = is small compared to the smallest time scale of the 
evolution problem. 

• All the phases, either solid or gaseous, have the same temperature held T{x, t). 
in] shows that this assumption it is valid even for the solid phase if the thermal 
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relaxation time tt^s = small compared to the smallest time scale of 

the evolution problem. 

Here we are interested in the mean behavior of a turbulent buoyant plume. Writing 
that solution we will use the following assumptions (see [m [ini DEI El [El [HI El El 

El ESI): 

• Reynold number is big enough and turbulence is fully developed, so that will 
be possible to disregard thermal conduction and shear dissipation. 

• Pressure is constant in horizontal section. 

• The profiles of mean vertical velocity and mean density in horizontal sections 
are of similar form at all heights. 

• The mean velocity held outside and near the plume is horizontal. We will need 
to make additional assumption on the dependence of the rate of entrainment 
at the edge of the plume to some characteristic velocity at that height. 

• Stationary how. 

• Radial symmetry around the source. 


2 The multiphase Dusty-Gas equations. 


Using the hypothesis given in the previous section, the Dusty-Gas model j2] 
simplihes: 


ApT + V • (Piu) = 0 , i e3 
+ V • (pju) = 0, j e a 
^tPm T V • {^PyPIL) 0 , 

+V-(^pmU®U-\-pi^= + Pm S' , 




■■ + Prau-g. 


(2.1a) 

(2.1b) 

(2.1c) 

(2.1d) 

(2.1e) 


As suggested in |T7] , it is convenient to use the specihc enthalpy hm = Cm + ^ = 
{Cra + Rra)T instead of the specihc energy Cm- We dehne the specihc heat at constant 
pressure of the mixture consequently: 


pi — (1 
^p,m 


+ Rn 


'^[yi{Ci + Ri)] + ) 

i& iea 


SO that hr 


— r T 


In this way, Eqs. (2.1) reduces to: 


V ■ {piu) = 0 , i G J 

V ■ {pju) = 0, j e a 

V ■ (pm w (g) w -f pi) = Pm S' 



( 2 . 2 ) 


(2.3a) 

(2.3b) 

(2.3c) 

(2.3d) 
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3 The Buoyant Plume Solution. 


Coherently with hypothesis of Section]^ we will look for a solution of Eqs. (2.3) 
in the following form: 


yk{r,z) 


Pmir.z) 


u{r, z) 


p(r, z) 
T(r,z) 


1, if r > b(z) and k = 1 

Ya(z) , if r < b{z) and k = 1 

0 , if r > b{z) and k ^ 1 

Yk{z) , if r < b{z) and k ^ 1 


(5{z ), if 0 < r < b{z) 
a{z ), if r > b{z) 


+U{z)z , 
-U,{z)f, 

< —Ue{r,z)f, 

Ue = Ue 

Mg —y 0 

V ^ 


if 0 < r < b{z) 
if r = b{z) 
if r > b{z) 
if r — b{z) 
if r 3> b{z) 


p{z) 

j Tfi{z), if 0 < r < b{z) 

|r„(z), itr> 6 (z) 


(3.1) 

(3.2) 

(3.3) 

(3.4) 

(3.5) 


where fc — i — 1 is the phase index corresponding to the atmospheric gas, while 
h ^4 1 is the generic index of a phase ejected by the plume vent. Here we used the 
so called purely “Top Hat” auto-similar profile. In general - as shown in PI - it is 

possible to use better prohles. Experiments show (see e.g. PI) that the auto-similar 
Gaussian prohle best ht data for a wide range of velocity measurements. Moreover, 
experiments are better reproduced choosing two different plume radius (say b{z) 
and Xb{z)) for the density and the velocity profile; the temperature profile should 
be determined by the equation of state of the fluid. Nevertheless, even if these 
modification could be done in Eqs. (3.1)-(3.5), here we decided - for simplicity - 
to use the “Top Hat” prohle. MISS(add comments to introduce Eq. (3.17) and 
the comments on entrainment, aggregation and settling that are included in the 
paper. Moreover, add something pointing out that we are neglecting the presence of 
humidity in the atmosphere) 

Here is an entrainment velocity. We shall write it as 




(3.6) 


where x is a dimensionless entrainment coefficient and is an arbitrary function 
of the density ratio (see e.g. 0). When ? 7 ^ = 1 we have the model of dH, if 
V>c{x) = \/x we get the model [H]. 

It useful to notice that inside the plume, the dusty gas constant and specihc 
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heat at constant volume Cm can be written: 


Rp = YaRa + = YaRa + > 


i=2 


i=2 


Cv,/3 — YaCy 


J:(Y,C.) + ZiYjCi) , 


(3.7) 

(3.8) 


i=2 


where Ra and Cv,q are respectively the gas constant and the specihc heat at constant 
volume for the atmosphere. We also dehne the specihc heat at constant pressure of 
the atmosphere and of the plume: 


Ca Cy^ci T Ra 1 

I 

Cf = + Rf = Y^c^ + (y,{Ci + Ri)) + E(yjC'j). 


(3.9) 

(3.10) 


i=2 


3.1 The mean conservation equations. 


For each altitude z E [0, L], we choose a control volume dehned as the cylinder 
of hxed radius B > b(z) centered above the source C = {{r, z) E [0, B] x [z, z + 6z]}. 


Using Eqs. (2.3a), (2.3b), (3.2) and (3.3), and the Gauss theorem, we hnd: 


0= / i:v.(ftt.) + yv.(p,u) = / v.(p„u) = 


, i& 


j&3 


= (3UttIP‘\z+ 5^ — l3Unb‘^\z — au^{B, z)2'kB5z . 


Now, dividing for 5z, sending it to 0 and then B -E b{z)^ we get total mass hux 
conservation: 


d.(0) = = 'iabU, . 


(3.11) 


In the general case, the source eject solid phases that are not in the atmosphere and 
some gaseous phase that is not included in the ambient composition. Identifying 
such a phases, respectively, with the index i E [2; /] and jGj = [/ + l;/ + J], and 


using again Eqs. (2.3a), (2.3b), (3.2) and (3.3), we hnd that the following mass 


huxes are conserved (we are neglecting particle aggregation and fallout): 

d.(g.) = d,{Y,l3Ub^) = 0 , Vz e [2; J], (3.12) 


d,iQ,)=d,iY,(3Ub^)=0, WjEd, 
while for the atmospheric phase i = 1 = a: 
d.(Q,) = d.(K„^W) = 2atf/.. 


(3.13) 


(3.14) 


Since the mass how rate of the erupted gases and particles are conserved, it is useful 
to dehne their mass how rate and mass fraction (respectively Qe.s and Ug^s): 


III I 

Qe = = Y.Qi = T.^ildub^ = qY.y, = qy,, 

i=2 i=2 i=2 i=2 

Q. = = Y.Qi = = QT,yi = Qy’ 

3 a 3 3 


(3.15) 

(3.16) 
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Putting together Eqs. (3.11), (3.12), (3.14) and 


/ i+j 

y^ + Y,y,+ Y. y, = y„ + Fe + h; = i, 

i=2 j=I+l 


(3.17) 


we obtain a relationship giving the mass flow rate Qa{Y) as a function of only vent 
conditions (Qi(0) = Qj(0) = Qyo) and Q{z)-. 


Qa{z) = Q{z) - I Y ) = Q{^) - (Qe - Qs) ■ (3.18) 

T=2 3 


By dividing Eq. (3.15), (3.16) and (3.18) by Q we obtain a relationship giving ns 
all the mass fraction as a function of only vent conditions {Qe,s) and the total mass 
flow rate: 


YJz) = 


Qe 


Q{z)^ 


YJz) = 1 - 


Qe T Qs 

~qW 


(3.19a) 

(3.19b) 

(3.19c) 


Dealing with the momentum, the vertical component of Eq. (2.3c) and Eqs. 


(3.2) (3.3) (3.4) yields: 


1 

6z 


-I3g = -TT^gJ - 7iag{B‘^ - J) = Y u^u + pz^ 


vr r. 


- - {l3U^J+pBQ 

dz L 


52^0 


Again, we take the limit B —)• b{z), obtaining 
dJ/3U^J) = {a-P)gJ. 


> dJnf3U b ) + tiB dzP. 

(3.20) 


(3.21) 


Here we used d^p = —ag, stated by Eq. (2.3c) together with p{r,z) = p{z) and 
n —)■ 0 when r b{z). 

Turning to the energy balance (2.3d) and using the same techniques, we And: 

(3.22) 




= 2abUf 


+ K\ -gfdUbJ 


where hp = CpTp and ha = CaTa- We neglect the term proportional to U^, to 
be compared to that proportional to t/^, because the entrainment velocity ?7e is 
typically one order of magnitnde smaller than U. 


Eq. (3.22) could be written in different ways using (3.11) and (3.21): 
d, {^UJCpTp) = {CaTj dzim") + ^dzim") - gaUJ , 


(3.23) 
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that is equivalent to Eq. (8) in [T7j, or 


d, {CpTp - cm) = d, (CM + - aoiUh^ , (3.24) 

where the dependence on the buoyancy flux and ambient stratiflcation is highlighted. 


Finally, we have that Eqs. (3.1)~(3.5) are one mean solution of (2.3) if 

'd^(Qe) = 0, 

dz(Qs) =0, 

d,(/3f/62) = 2abU, (3.25) 

d ,(/ 3 f / 262 ) = ( a -/?)^&2 

[d, - CM) = d, {CM + ^dMb^) - gc^ub^ • 


By noting again that Qe and Qs are conserved and that Eqs. (3.19) hold, here the 


unknowns are /3{z), U{z), b{z) and Tp{z)^ provided the knowledge of the ambient 
density a, the ambient temperature Tq, and the dependence of 17e on the other 
unknowns (the entrainment model). We are still lacking in one condition. The 
equation of state of the various phases together with the full expanded plume 
hypothesis - p{r, z) = p{z) - will give us that last needed condition. 


4 The Gas-Particle Plume model. 


In order to close the latter system of equations, we can use solution (3.1)-(3.5) 


with the constitutive law for the dusty gas pressure. Since in Eq. (3.4) we have 


assumed p(r, z) = p{z), we have that - at a given height - the pressure inside the 
plume is the same of that outside the plume: 

p = (3RpTfi = aRCT^. (4.1) 

Thus, we can rewrite the plume internal-external enthalpy differential as follows: 

RaCfS 


MpTp - CM = aCMif^ - MTa . 

KpUa 


(4.2) 


We define the thermodynamic properties of the ejected gas and of the particles as 
follows 

(4.3) 

e i=2 
1 ^ 


Rc = Y ^ 


Ce=Y Y."d^{C^ + R^), 


e i=2 


Cs = Yj:y^c,, 


(4.4) 

(4.5) 


noticing that all these quantities are - coherently - conserved along In this 
way thermodynamic properties of the mixture can be written in terms of the 
thermodynamic properties of the three components, for example: 

Cf, = YM + Y,C, + Y,C,. (4.6) 


^It is sufficient to multiply both numerator and denominator of the right hand sides by Q, 
and notice that YkQ = Qk = Qk,o- 
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Using these definitions pins Ys = Ye = Ue = p , and Eqs. (3.17), (3.15), 

Oa -^ O ! I I ' I I ' 

(3.16), we can write in a convenient form Eq.(|4.2|): 

- CM = ccr, 


{a — j3) + a 

Now, defining the relative flux of enthalpy 

X.sQs (Xe '4’e}Qe 


XsQs T (Xe '4^e)Qe 


F = 


{a — (3) + a 


{Q — Qs) + (V'e ~ l)Qe 


(Q — Qs) + (V’e ~ l)Qe 


[/F 


(4.7) 


(4.8) 


equation (3.24) can be rearranged 


CnTn 


2CaTa CaTa 


(4.9) 


It is useful to define 


Qip Qs 3~ ("06 I)f^e : 

Qx ~ ~ l)Qs + (Xe ~ l)Qe ) 


which are constants along z, so that 


F 



(3) F oi 


Qx Qip 

Q T Qpj 


Ub"^. 


(4.10) 

(4.11) 


(4.12) 


This expression for F represents a modification of the buoyancy flux for a dusty- 
gas plume in the general non-Boussinesq case (cf. [Hj). It takes the classic form 
(a — (0.0) for a single-component gas plume (in such a case Q^ = 0 and 

Qpj = 0). For this reason we will refer to the relative flux of enthalpy F as the 
dusty gas buoyancy flux, a generalization for the multiphase case of the standard 
buoyancy flux. 

This new quantity F, together with the mass fliix Q = /3Ub‘^ and the momentum 
flux M = (3U%‘^ allow us to close problem (3.25) in their terms: 




Q’ = 2U,(a,Q,M,F). 

M [ E 

OL 

aQ{F + Q){Q + Qip) 

M[Q + Qx] 

Q){Qx ~ Qb)l 
"[Q + Qx] \ 

' M^Q' g{F + Q){Q + Q^) 

2C'„T„Q2 CMQ + Qx) ’ 


(4.13a) 

(4.13b) 

(4.13c) 


where U = ^, b 


Q{F-\-Q)(Q-\-Q^) 


and fl = a 


QIQ+Qxl 

{F+Q)(Q+QF 
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5 Non-dimensionalization. 


It is useful to transform the latter problem in dimensionless form. We choose 
Q{z) = Qog(C), M{z) = Mom{C), F{z) = Fof{C) and ^ = 4C (4 = where 

(■)o refers to the vent height. In this way, we have g(0) = m(0) = /(O) = 1. It 
is worth noting that C = 0 can correspond to the actual vent elevation as to any 
height above the vent (cf. [5]). The model in non-dimensional form then is 


q = VgT]^ 


m{(pf + q)iq + q^) 


m = Vm — 
m 



f = 


taiO 


qiq + qx) 
i<Pf + q) \ 
{q + qx)) 


(0/ + g) 0/(0 - 


q+jiA (t> rn^q' 

q + qx) q 


(5.1a) 

(5.1b) 

(5.1c) 


where - dehned in Eq. (3.6) - is the entrainment function, potentially depending 
on the other variables and parameters; a{() = a{£oC)/ao, ta{C) = Ta{ioC)/Ta,o, 
0 = Fo/Qo, = Q^/Qo, qx = Qx/Qo, 7c = 0/(C) = -^^a(C) and 


Vq = 2 X 


'^m. — 


= 


gFoQoio ^ 
■ 


Mi 

gQo£o 


= Ri 


gh 


g^o 


gio 




(5.2) 

(5.3) 

(5.4) 


Fq CaTafl (j)CaTafl CiSflTpfl — CaTafl A/Iq 

We call these last three parameters the rate of variation respectively of q, m, f. 


In Eq. (5.3), we have given a modihed dehnition of the Richardson number Ri = 


(pgio/Ui, because cjyg = g' in the monophase case {g' being the reduced gravity). In 


Eq. (5.4) we used the dehnition of the Eroude number Er = Ui/gio and of the Eckert 
number Ui/Aho, where Aho = CpflTp^Q — CaTafl is the enthalpy anomaly at the vent. 


Moreover, we have used Eqs. (4.7), (4.8) implying (j)CaTafi = CpfiTjsfi — CaTafl . It 
is also useful to rewrite the physical variables as a function of these new parameters: 


Mq m 

Qo q 


b 


^ q{4>f + q){q + qj^) 

am{q + q^) 


/3 = a 


qjq + qx) 

(0/ + O(9 + 0y) 


T0 = 

h)j(s) 


rj. ^f + q 

J-OL 

q + qx 

h)5,o(s,o) 

q 


(5.5a) 


(5.5b) 


(5.5c) 

(5.5d) 

(5.5e) 


It is worth noting that q^ > —1 because the specihc heats and gas constants 
are positive (y., ip. > 0) and the sum of the initial mass fraction is smaller than 
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1 (cf. definition of in Tab. 5.1). Moreover, 0 > —1 becanse > 0. 


Even if these are the general conditions for such parameters, in Tab. |5.1| there are 
summarized the possible ranges for volcanic eruptions. 

Using d^p = —ag and the ideal gas law it is possible to obtain the density 
stratification as a function of the temperature: 


a(C) = ta (C) exp 


9^0 


K 


RryT, 


a -^ Q :,0 *^0 


t-AC)AC' 


(5.6) 


For example, if the non-dimensional atmospheric thermal gradient 6 = Oaio/Tafl is 
constant, we have ta{C) = ^ ~ 


a{() = {i-ec)^< 


--1 


(5.7) 


and 0f{() = Of = O/vfcj). 

It is also useful to define the Brunt-Vaisalla frequency N. Recalling that the 
potential temperature is 


t^A0 = U0 KC)ta(C))“^, 


we obtain 


= l ln(ip,„)'(C) 

% 


9^ l-OfiO 

taiC) 


CaTafi 


(5.8) 


(5.9) 


This frequency depends on the height z, but it can be approximately be considered 
as a constant because it vary slowly in our atmosphere: ~ 10 % of variation in the 
troposphere. In what follows we call Nq its constant approximation. Using standard 
average conditions for the troposphere, we find Nq ~ 1.13 * 10“^ Hz. Studying 


plumes in a stratified atmosphere (cf. Sec. 5.6), it is useful to define 


l-9f M 






to 


^9 ^9 


= Vffi , 


(5.10) 


showing that the new parameter Vf^ can be recovered by knowing the enthalpy 
anomaly 0 and the non-dimensional stratification length scale ^ = g/NAo- In other 
words, the more Vf^ increases the more the vent dimensions corrected with the 
enthalpy anomaly are comparable with the stratification length scale. 

All these non-dimensional parameters characterize the multiphase plume and give 
us the possibility to classify through them all the possible regimes. We summarize in 
Tab. ^ six of them, which are the independent non-dimensional parameters 
snfRcient to characterize a multiphase plume. In order to fix ideas, we show 
there the range of variability of those independent parameters for Strombolian to 
Plinian volcanic eruptions. 

Indeed, the knowledge of these parameters and of the thermodynamic properties 
of the atmosphere allows us to retrieve the physical dimensional parameters. We 
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parameter 

explicit form 

range of variability 

description 

</> 

C/sTiSfi — CaTafi 
CaTafl 

0.3 5 

enthalpy anomaly 
(non-Boussinesqness) 

gp 

~+s,0 + (06 ~ l)le,0 

-1 1 

mass flux anomaly 
due to gas constants 

(lx 

.Xs — l)+s,0 + (Xe — 1)1^,0 

-1^1 

mass flux anomaly 
due to specihc heats 

vj2 

K 

0.05^0.3 

entrainment 

coefficient 


(pgio 

ui 

10“^ 10 

modified 

Richardson number 

i 

9 

Niio 

10^ ^ 10^ 

stratification 

length-scale 

Tab. 

5.1: Independent parameters for a 

multiphase plume in a stratified atmosphere. 

report here all the inversion relationships needed: 


if) = S-^. see footnote 

Nii 


(5.11a) 

O 

o 

(1 + 0)(1 + q^) 

1 + 


(5.11b) 

o 

o 

1 + 


(5.11c) 

(1 + 0)(1 + q^) 


V Vm 


(5.11d) 

rr rr 1+0 

Tm - r,,„ J ^ 


(5.11e) 

Qo = f^oUobo 


(5.11f) 

Mo = f^oUX 


(5.11g) 

Fq = 0Qo 


(5.11h) 

(lx 

7c = 

- 9^ 

0 


(5.11i) 


9^0 


(5.11j) 

0C'„T„,o 



Vf,o = 
Ye,0 = 


Qx + (Xs - i)g^ 


(Xe - 1) + (Xs - l)(V’e - 1) 

^,0 (^e 1)^,0 Q'lp 

Ya ,0 = 1 ~ ^ s ,0 ~ ^,0 • 


see footnote 


(5.11k) 

(5.111) 

(5.11m) 

(5.11n) 
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In Cerminara et al. [3], we have used these inversion relationships to obtain the 
vent condition of a real volcanic eruption occurred at Santiaguito (Santa Maria 
Volcano, Guatemala). 

In this thesis, we will study only two of all the possible entrainment models 
introduced in the literature: 


Morton et al. m, where = 1 


Ricou and Spalding where = rj^[(3/Q.) = 


QiQ + Qx) 




More complex models have been studied in volcanology and fluid dynamics. One 
example can be found in pQ where depends on the local Richardson number. 

It is worth noting that the mass flux q{() is a strictly increasing function as 
long as 7]k is positive, while the sign of m'{() depends on the buoyancy sign: 


sign (buoyancy) = sign / — 7< 


( 0 / + ?) ' 

+ (lx), 


(5.12) 


because Vm, Q, "m are strictly positive. For an analysis on the plume buoyancy 


behavior see Sec. |5.4[ In Sec. |5.6| we will study in detail the evolution of the plume 
variables under the Boussinesq approximation. However, something can be noted 


even at this point of the analysis by looking at the full system (5.1): 1) the mass 


flow q(z) is a strictly increasing function because the entrainment models we are 
using are positive functions; 2) the momentum flux m{z) has derivative equal to zero 
when the buoyancy become zero. It can be due to two causes, buoyancy reversal 
or neutral buoyancy level. We denote Cnbi the neutral buoyancy level; 3) when 


m{z) = 0 system (5.1) encounters a singularity. In that point the plume reaches 
its maximum height Cmaxi 3) the enthalpy flux is a strictly decreasing function, 
because usually in applications the term containing is dominant and 

negative. 

In the next sections we discuss some of the approximations applicable to problem 


(5.1). In particular we find that 7c is the parameter related to the column instability 
- if 7 c > 1 then the volcanic column will collapse - and that (j) is the parameter 
measuring the non-Boussinesqness of the mixture - if 0 1 then the Boussinesq 

approximation holds. Moreover, q^ and q^ are the parameters measuring the 
multiphaseness of the mixture - if |g^| ~ -C 1 the plume can be considered as 
a single phase one. 

In this thesis we will study three different volcanic eruption and one experimental 
plume that we denote, from the weaker to the stronger: [forcedPlume], [Santiaguito], 


[weakPlume], [strongPlume]. We report in Tab. 5.2 all the parameters for these 
volcanic eruptions, respectively: 1) the physical parameters at the vent - radius, 
density, temperature, velocity and mass fractions; 2) the mass, momentum and 
enthalpy flows; the non-dimensionalization length scale and the multiphase Morton 


^When stratification is disregarded, no reference length scales are present in the non- 


dimensional system, thus 6o must be given and fo can be recovered from Eq. (5.11b I. 

^In order to have the mass fraction of ejected gas and solids, their thermodynamic properties 
must be known: namely their specific heat and the gas constant of the ejected gas. 
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length scale (see below); 3) the six independent non-dimensional parameters; 4) the 
non-dimensional dependent parameters; 5) the non-dimensional plume maximum 
and neutral buoyancy level height, as obtained from system (5.1) with Ricou and 
Spalding [H] entrainment model 


5.1 Monophase plume. 

If the thermodynamic properties of the ejected fluid are similar to them of the 
ambient fluid then ~ \q^\ -C 1. In this case, model 


(5.13a) 
(5.13b) 
(5.13c) 

where 


q' = Vgr]^^a{z) 

qf 


+ q) 


m = Vr, 


m 


f = 


_V_ 

ta{z) 


W+ q)iOfiz) - 1 ) -f- 


(j) m^q' 
2nm q^ 


(5.1) becomes: 


= 1 


(im) 

m)- 


It is worth noting that in the single phase case = Ca and Rp = Ra- Thus, 
the initial enthalpy anomaly reduces to the initial thermal anomaly or equivalently 
to the density anomaly: 

, _ T)3,o ~ Taft _ ATo _ ao — /^o 

Ta,0 Tafl /^O 

Consequently the reduced gravity becomes g' = (pg. 


(5.14) 


5.2 Jet regime 


In the jet regime - defined as the one where m = / = 1 - Woods na pointed 
out that the Ricou and Spalding [T3] model can be used. In this case, Eqs. (5.1) 
simplify a lot, becoming: 


q 


f 


= Vn 


m' = 0 


/' = 0, (5.15) 


with the easy solution q{() = VgC + 1. 

Substituting this solution in Eqs. (5.1) and proceeding with the dimensional 
analysis, it is possible to hnd Im, the dimensionless transition length scale between 


■^While for [forcedPlume], [Santiaguito], [weakPlume] we have used a constant atmospheric 
thermal gradient, for [strongPlume] the atmospheric temperature profile is a little bit more complex, 
because we have included in it the presence of the tropopause [cf. 0]. 
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parameter 

[forced Plume] 

[Santiaguito] 

[weakPlume] 

[strongPlume] 

bo [m] 

0.03175 

22.9 

26.9 

703 

f3o [kg/m^] 

0.622 

1.05 

4.87 

3.51 

ao [kg/m^] 

1.177 

0.972 

1.100 

1.011 

Tp,o [K] 

568 

375 

1273 

1053 

T«,o [K] 

300 

288 

270.92 

294.66 

Uo [m/s] 

0.881 

7.29 

135 

275 

[mVs^K] 

287 

287 

287 

287 

[mVs^K] 

1004.5 

998 

1004 

1004 


- 

1.61 

1.61 

1.61 

Xe 

- 

1.866 

1.803 

1.803 

Xs 

- 

1.102 

1.096 

1.096 


0 

0.196 

0.03 

0.05 


0 

0.410 

0.97 

0.95 

ya,o 

1 

0.394 

0 

0 

No [Hz] 

1.14* 10-2 

1.43 * 10-2 

1.40* 10-2 

2.14* 10-2 

7tQo [kg/s] 

1.74 * 10-3 

1.26* 10^ 

1.5 * 10® 

1.5 * 10*^ 

ttMq [kg m/s^] 

1.53 * 10-3 

9.19* 10^ 

2.02 * 10® 

4.12* 1033 

vrFo [kg/s] 

1.55 * 10-3 

7.28 * 103 

6.35 * 10® 

4.56* 10^ 

io [m] 

0.02308 

23.8 

56.6 

1310 

[m] 

0.0854 

18.4 

352 

4070 

0 

0.893 

0.58 

4.25 

3.04 

gp 

0 

-0.290 

-0.952 

-0.920 

<lx 

0 

0.212 

0.117 

0.131 

Vq 

0.28 

0.659 

0.2 

0.2 

'^m 

0.261 

2.54 

0.129 

0.517 

1 

3.29 * 10® 

2020 

886 

16.4 

7c 

0 

0.869 

0.252 

0.345 


8.41 * 10-^ 

1.41 * 10-3 

4.81 * 10-^ 

1.43* 10-2 

^f,0 

3.40 * 10-^ 

8.58 * 10-^ 

2.66* 10-^ 

2.01 * 10-2 

Cmax 

1665 

23.98 

160.6 

20.68 

Cmax/ Cnbl 

1.318 

1.306 

1.354 

1.523 


Tab. 5.2: Relevant parameters of the plumes studied in this thesis. 
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the jet and the plume regime. It is the length scale for which the momentum 
variation becomes important. From the momentum equation we hnd: 


1 


— ~ VmiVqiu + 1 ) ^ VmVqi 


M 




Km = [VqVm) 


(5.16) 


from which, back to dimensional units: 

r ( UUO \ ^ 


(5.17) 


This quantity became equivalent to that defined in Morton HDl when = qx = ^ 
and /3 ~ a. 

The typical length scale of stratihcation for a jet can be found by using a 


similar dimensional analysis for Eq. (5.1c) 


or 


IT = 1 + Vs) - 'VfflVq^s ^ 4 = (w,o) " , 

ts 


(5.18) 


“ U,f,oy ~ UoN~ 


(5.19) 


This parameter is comparing the rate of variation of m and /. We have that if (5j < 1 
than stratihcation have a role in the jet-like part of the plume, on the contrary, 
if > 1 stratihcation is important just in the plume-like part of the plume. We 
will comment better this length scale in the section below dedicated to the plume 
height. 

Usually in jets, atmospheric stratihcation is not important because of their 
limited height ((5j > 1). We want to explore now when the kinetic correction term 
could be important. Contrarily to the last two terms, the second term in square 


brackets in Eq. (5.1c) become less important as ( grows. In particular it decreases 
with q'/q"^ cx Dehning the typical length scale for this term we have: 


ix 2Vm{l+Vq£Ky 





(5.20) 


admitting a positive solution if and only if 


Vm, 


( 4:1) 

Ur) 


(-K ^ < 


Vq(pVf 

(t>Vf 

‘2VqVm 


2Aho 

IP 

AxAHq 


1 

> 1 . 


(5.21) 


Thus, the kinetic correction can be important just near the vent or very far from 
it and only when Aho -C Uq (Ec S> 1). In other words, this correction can be 
important for “cold and fast” jets and far from the jet center. Generally, in volcanic 
plumes the Ec number is small, thus the kinetic correction can be disregarded. 
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5.3 Non stratified plume regime 


If stratification and the last term in sqnare brackets of Eq. (5.1c) can be 
disregarded, / = 1 and model (|5.1|) becomes 


q = VgT]^ 


/ ^ / 1 
m = Vm— 1 


m(0 + g)(g + g^) 

^| qi.q + qx) 

(0 + g) \ 


m 


(g + gx). 


/' = o 


(5.22a) 

(5.22b) 

(5.22c) 


This ordinary differential eqnation has a hrst integral of motioi:[^ IX in both the 
considered cases for rj^. We fonnd respectively for the entrainment models of 
Morton et al. m and Ricon and Spalding El- 


u 


MTT 


= 2 




(0 + g) \ g(g + qx) 

(g + gx)/ \ (<^ + g)(g + gy) 


dq — 


Rrs = g^(l -7c) - 27c(0- gx) [g - gxln(|g + gxl)] 


5nm 

4y 


5n. 




(5.23a) 


(5.23b) 


Using this hrst integral of motion in Eq. (5.22a) it is possible to hnd an implicit 
solntion for the height of the form C, = C{q)- For the Ricon entrainment model, 
dehning 

/(g) = g^(l - 7c) - 27 c (0 - gx) [g - qx ln(|g + gxl)] , (5.24) 


and snbstitnting the corresponding hrst integral of motion fonnd in Eq. (5.23b) 

dv. 


URs(g,m) = /(g) 


5r;r 


-m5/2^URs(l,l)=/(l) 


4ng 

5Vrr 


into Eq. (|5.22a|), we fonnd the following implicit solntion: 

5n^ 


C = C(g) = 


'1 


7.1 

V, 


Q L 


4:Vn 


(/(x) — Rrs) 


(5.25) 


(5.26) 


Using this solution it is possible to hnd the height at which the Boussinesq approxima¬ 
tion starts to hold: ( = Cbou- We choose the value q = gsou = lOmax(|0|, |g^|, |g^|). 
In Tab. 5.3| are reported the value we obtain for the examples considered in this 
thesis. By comparing those values with (Cmax reported in Tab. 5^ it is possible to 
have an idea of the part of the plume where the Boussinesq regime holds. 

Under the same hypothesis of this section, the monophase case (5.13) becomes 
equivalent to the model studied in Fannelpp and Webber [5]: 


g ='Vqr]^i 


+ q) 


m = Vr 


f = 0. 


<1 

m 


(5.27) 

(5.28) 

(5.29) 


first integral of motion is a quantity remaining constant along the motion described by the 
differential equation. It is also called constant of motion. 
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For the entrainment models of Morton et al. m and Ricou and Spalding m the 
hrst integral of motion are respectively: 

Umtt = (? - \/?(9 + 0) + In + ^Jq + (i^ - (5.30) 

Rrs = . (5.31) 


5.4 Buoyancy reversal and plume stability 


In this section, we consider the plume model behavior near the vent, where it 
is not possible to use the approximation g S> | 0 |, Ig^l) |gpl (see next section) but 
/ ~ 1 as done in the previous section. Here we will use the Richou entrainment 
model because we are near the vent, however the present analysis is independent 
from the entrainment model used since the sign of the buoyancy does not depend 
on r^K- In model (5.22), the sign of the buoyancy force is determined by: 


sign (buoyancy) = sign 1 — 


(0 + g) ' 

(g + gx). 


= sign(/'(g)) 


(5.32) 


Here, l{q) is the hrst integral function dehned in Eq. (5.24). When l'{q) < 0, 
the plume is negatively buoyant and m decreases. If we arrive to the condition 
l{q) = Rrs then m —)■ 0 because the hrst integral Urs niust be constant. Thus the 
plume stops (or collapses) and it is not able to reverse its buoyancy. 

We can better understand the behavior of the non-stratihed multiphase plume 
by analyzing all the possible conhgurations. For this purpose, it is useful to dehne 


* ^ 1 + gx 
1 + 0 


(5.33) 


To,a _ bc0 gx 

7 ti ) gmin ) 

+),/? 1 - 7c 

where /^gmin) = 0. We enumerate the following situations for g > 1 (recall that 
g(C) ^ 1 because it is a strictly increasing function and 0 > — 1) by denoting 
“C” the cases when the plume collapses and by “B” the cases when the plume can 
reach and sustain the condition of positive buoyancy: 

IB) positive buoyant. If 7 c < 1 A 7 c < 7 * 

then /'(g) >0 Vg > 1 and the plume rises indehnitely. 


2B) zero, then immediately positive buoyant. If 7c = 7 * < 1 
then/'(g) > 0 Vg > 1, /'(g) = 0 if g = l 


(t>>Qx) 


3B) jet with zero buoyancy. If 7 c = 7 * = 1 (^ 0 = g^) 

then /'(g) = 0 and the plume behaves as a jet. 

4BC) from negative to positive buoyancy. If 7 * < 7 c < 1 (=^ 0 > q^) 

then /'(g) < 0 when g < gmin, the minimum of /(g) is reached in g = g ^in 
and /'(g) > 0 when g > gmin- In this case inversion of the buoyancy 
sign can be possible if the minimum value of /(g) is above the hrst integral: 
/(gmin) — Rrs > 0. In the opposite situation /(gmin) — Rrs < 0 the plume 
is not able to invert its buoyancy and it collapses when m = 0 , thus when 

/(g) = Urs. 
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parameter 

[forcedPlume] 

[Santiaguito] 

[weakPlume] 

[strongPlume] 

Cbou 

15.77 

7.07 

82.2 

53.9 

7c 

0 

0.869 

0.252 

0.345 

7* 

0.528 

0.768 

0.213 

0.280 

Qmin 

- 

2.22 

1.27 

1.40 

/(<?min) UrS 

- 

0.0388 

1.19 

0.214 

Oq 

0.860 

1.59 

1.65 

0.473 


Tab. 5.3: Column stability parameters for the plumes studied in this thesis. 


5C) from positive to negative buoyancy. If 1 < 7c < 7* 0 < (lx) 

then l\q) > 0 when q < q^i^n the maximnm of l{q) is reached in g = g^in 
and /'(g) < 0 when g > gmin. In this case the plume always collapses going 
from positive to negative buoyancy. 

6 C) zero, then immediately negative buoyant. If 7 c = 7 * > 1 0 < g^) 

then/'(g) < 0 Vg > 1, /'(g) = 0 if g = l 

7C) negative buoyant. If 7 c > 1 A 7 c > 7 * 

then /'(g) <0 Vg > 1 and the plume collapses being always negative 
buoyant. 


Thus, we can summarize that: 1) if 7 c > 1 the plume starts or becomes negative 
buoyant and collapses; 2 ) 7 * must be compared with 7 c to know the initial buoyancy 
of the plume: if 7 c < 7 *(>) then the plume is initially positive (negative) buoyant; 
3) if 7 c < 1 then the plume is or can become positive buoyant, buoyancy reversal 
occurs if /(gmin) — hiRs > 0. In Tab. 5.3 we report all of these parameters for the 


plumes studied in this thesis. While [forcedPlume] is positive buoyant, the other 
three plumes are initially negative buoyant. For all of them, buoyancy reversal 


occurs. 


5.5 Non stratified Boussinesq regime 


In the Boussinesq limit, we have that g 3> |0| , |gxl ! It is worth noting that 
under this approximation the reduced gravity g' can be written via (/>: 


(fg ~ 


tto — /3o 

tto 


9 = 9 - 


(5.34) 


Moreover, the two entrainment models we are considering become equivalent and 
Eqs. (5.1) reduces to: 


g' = Vq\/rn 
m = Vm (1 
/' = 0 . 


7c) 


m 


(5.35a) 

(5.35b) 

(5.35c) 
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which is the multiphase version of the celebrated model introduced by Morton et al. 

m- 



/ q 

m = Vm — 

m 

f = 0. 


(5.36a) 

(5.36b) 

(5.36c) 


Thus, we have found that the equations for a multiphase plume in a calm environ¬ 
ment under the Boussinesq approximation are equivalent to the monophase Morton 
et al. HH model with the following modification: 


Vm -t Vmil -7c) • 

Model ( 5.35[ ) has the following first integral: 


Umtt — Idns — U — q 


^Vm{l - 7c) 


m 


5/2 


It = 1 — a„ 


AVn 


dq — 


5^;„^(l - 7c) ’ 


(5.37) 


(5.38) 

(5.39) 

(5.40) 


The values of for the plume examples studied in this thesis are reported in 
Tab. 5.3 From this expression and Eq.( 5.35a), we found the implicit solution: 

C = C{q) = [ dx x‘^ — 1 + ttq ® . (5.41) 


This solution has two branches, depending on the sign of (1 — 7c), thus on the sign 
of ttq. If ttq < 0, the column is unstable with implicit solution (cf. App. [^for the 
definition of the Gaussian hypergeometric functions '^b and 0?,): 


c = 




Vq{l - aq)5 






(5.42) 


The maximum height is reached when q^ax = \/l — a^: 


Hmax/ 


(-ag)i 


Vq{l - aq)5 




1 - a„ 


(5.43) 


In Fig. 5.1 we show the behavior of i7max/f*o foi' Vq = 0.2 and we compare it with 
the following asymptotic expansion (t5^-i/5(1) — 1.150): 


Hraax /^0 — — 1^+0 (( —Ug) 


(5.44) 


Thus, the maximum height of a collapsing multiphase plume in Boussinesq regime 
behaves approximately as y/—d^- 
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Fig. 5.1 : The height of collapse of a multiphase plume in a non-stratified stable atmosphere as a function of the 
parameter aq defined in Eq. | |5.40[ l. Here we compare the exact formula Eq. | |5.43[ | with its asymptotic 
expansion Eq. | |5.44| l, in the case Vq = 0.2. 


On the other hand, if aq > 0, the column is stable, rising indefinitely with this 
law (see App. 0: 





(5.45) 


The asymptotic expansion (5{x) = l + 0{x) allows us to find the self-similar solution: 


?(C) = 



m(C) 


LLq 


{q\C) - i) +1 


oc 


(5.46a) 

(5.46b) 


From here it is possible to extract the asymptotic plume radius evolution: 

^(0 =-1= = |^^gC + a|<S_i(l-s). (5.47) 

y^m(C) ^ 

In this formula, we can recognize the famous result of Morton et ah m- the plume 
spread b'{Q is asymptotically constant and equal to ^Vq = |x- Moreover we found 
the initial virtual radius of the asymptotic plume and its asymptotic approximation. 


K = 1 (1 - aq) ~ 0.5012^^;) + 0.6 . (5.48) 

The virtual plume radius is the intercept between z = 0 and the radius of the 

g^l/5 

equivalent plume spreading from a point source at z = Zy = —~ 
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(a) Virtual radius 6v (b) Necking height Cneck 

Fig. 5.2 : a) The virtual radius bv as a function of aq. The virtual radius tends to zero when aq 0 and increases 
with a square root law as aq increases (cf. Eq. ( |5.48[ l). b) Height of the plume radius necking Cneck 
predicted by Eq. l |5.51| l. 


In Fig. 5.2a we show the behavior of by{aq) and of its asymptotic approximation. 


Finally, it is worth noting that the derivative of the plume radius has a simple 


expression thanks to the first integral (5.38) 


b'iC) = Vq 


from which 


2(1 


baqUi^/'^ 


(5.49) 


b'{0) = Vq 


3 2(1-g,) 

5 5a„ 


(5.50) 


is the plume radius slope at C = 0. Another important property is the necking 
height C = Cneck, where &'(Cneck) = 0. It exists only when 0 < < 2/5: 


Cneck — „ 

3v„ 


3(l-«g) 


<S. 


-) - 0_i (1-a,) 


(5.51) 


As shown in Fig. 5.2b the necking height never exceeds C = 1- 


We summarize in Fig. 5.3 all the possible regimes of model (5.35). Ranging from 
Qq = 0“ to ttq = O’*" passing through Qq = oo, we have shown that: 1) {collapsing 
regime) when Oq < 0 the plume is collapsing, b'(0) > Vq, and its height increases 


as Gq decreases (cf. Fig. 5.1); 2) (jet regime) when Og —)■ cx) then model Eq. (5.35) 


reduces to the jet model (5.15) with b{z) = VqZ + 1; 3) (forced plume regime) when 


> 1 the initial slope is ^ < 6'(0) < Ug, and the plume starts behaving as a jet 


until z < I'm (cf. (5.16) and Morton (lU]), then it moves to the plume-like behavior. 
As shown in Figs. 5.3 5.2a, and b^ increase with Og; 4) (pure plume regime) when 
Og = 1 the solution of model (5.35) highly simplihes and asymptotic expansions 
coincide with the exact solution. In particular, we have b(z) = ^z + 1. There is 
not a jet-like interval in this regime; 5) (buoyant plume regime) when 0 < Og < 1 
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. collapsing: < 0 -jet: a^ = ^ -forced plume: l<a^<oo 

— ■ — pure plume: = 1-buoyant plume: 0< a^<l 


Fig. 5.3: Evolution of the plume radius h{z) = qjy/in in all the admissible regimes of model l |5.35[ | with Vq = 0.2. 
Starting from the lower graph, we choose: aq = —1, —10, —50, co, 50, 10, 1, 0.1, 0.0001. 


we have 6'(0) < and the plume radius reach its asymptotic slope ^ rapidly, 
after a small necking interval. In particular, if 0 < < 2/5 there exist (/neck > 0 

where 6'(Cneck) = 0. 


5.6 Boussinesq plume regime in a stratified environment 


The Boussinesq approximation, with atmospheric stratihcation reduces (5.1) to: 
q' = Vq^Ja{C,)m 


m' = Vm-U-lc) 


m 




(5.52) 

(5.53) 

(5.54) 


If we consider the atmospheric stratihcation only at the hrst order, we can apply 


the following approximation to the latter system (cf. Eqs. (5.9) and (5.10)): 


a(C) ^ 1 




1 - OAO 
taiC) 


- ^/,0, 


(5.55) 


allowing us to write the multiphase plume model in a stratihed calm atmosphere: 


q = Vq\/m 


m' = Vm — if - 7c) 


m 


f = -Vf,o q ■ 


(5.56a) 

(5.56b) 

(5.56c) 
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This model reduces to the same model introduced by Morton HD] in the monophase 
case: 


q = Vq^/m 

(5.57a) 

/ Qf 

m = Vm — 
m 

(5.57b) 

' O 

1 

(5.57c) 


where Vf^ is proportional to the Brunt-Vaisalla frequency iV^ (cf. Woods [IH] and 
Eq (|5T0|). 


In order to find the first integrals of motion, we write system (5.56) in this form: 


dq m dm df 

Vqy/rh Vmqif - 7c) ’ 


(5.58) 


By using the last equation multiplied by q{f — 7 c), we obtain the first conserved 
quantity (recall that /o = mo = 1): 


U 


m 


'Vm 


+ {f - IcY 


(l-7c)^ + 


'^m 


(5.59) 


Urn is a very interesting quantity, because it holds whatever the entrainment model 
is. Indeed, we have found it just by using the conservation of mass and enthalpy in 
system (5.56), which are independent from the entrainment model. Moreover, this 
conserved quantity is telling us that m reaches its maximum value 


'^max 



V-m 

^/,0 


(1 - 7c)2 , 


(5.60) 


when / = 7 c. In other words, the flux of momentum is maximum when the flux of 
buoyancy (/ — 7 c) is zero: neutral buoyancy level. 

Additionally, this first integral of motion tells us the value of the enthalpy flux 
when the plume reaches its maximum height. We define the maximum height of the 
plume as the point ( = Cmax where m = 0, thus the minimum value of the enthalpy 
flux should be 


/ (Cmax) — /min 7c 


(5.61) 


because / is a strictly decreasing function of ( (cf. Eq. (5.56c)). Thus, increasing 
the height ( from 0 to Cmax let / decrease from 1 to /min; while m increases from 1 
(/ = 1) to rumax (/ = 7c), then it decreases to 0 when / = /min- These observations, 
will be very useful in the next sections of this chapter. 


Moving back to Eq. (|5.58|), it is easy to show that: 
qdq = 


-\/mdf = 


V 

UqUm 


(«.»-(/-70)")'" d/, 


5/4 


Wo 

from which we obtain another first integral of motion: 

/(/-7c)'\ 


(5.62) 


Uq — q -|- 


5/4 

-^/.O 


U/‘‘(/-7c)5! 


u. 


(5.63) 
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parameter 

[forced Plume] 

[Santiaguito] 

[weakPlume] 

[strongPlume] 

1 -Xo 

6.521 * 10“^ 

9.163* 10“® 

1.832 * 10“® 

4.245 * 10-2 

% 

2.054* 10“® 

1.753* 10-® 

2.178 * 10“^ 

3.940 * 10-2 


1.114* 10“® 

0.1363 

6.062 * 10-2 

0.3010 

Op 

0.9321 

0.5183 

0.4828 

1.691 

Cmax 

1526 

21.35 

141.9 

17.89 

x(i) 

Smax 

1524 

20.54 

139.3 

15.89 

m 

Smax 

1532 

24.16 

151.5 

24.33 

Cmax/ Cnbl 

1.318 

1.375 

1.345 

1.488 

di/dbl 

1.318 

1.394 

1.354 

1.582 


Tab. 5.4: The main parameters defined in this section for the four plume examples of this thesis. 


where 5^i(a;) = 2 F 1 is the hypergeometric function dehned when x <1 

in App. B and ~ 0.874C®, Noting that is a 

4 _ 4 

strictly increasing function bounded in [—1,1], we have that, as / decrease from 1 
to 7 c — ^/Um, q must increase from 1 to 


?max - 1 + 




■'q^m -If 1/4 
5/4 
^/,0 


(l-7cm 


' (l-7c)^ 

If™. 


+ V him 5^1 (1) 


(5.64) 


problem (5.56): 


By using again Eq. (|5.56c|) with (|5.63|), we have found the implicit solution of 

'if-lor 


C = — df 

Vf,oJ 


ll„ 




q^m 

5/4 

A,o 


-Ki‘ (/' - 7c) 


'Ll 


(5.65) 


In order to better understand the behavior of the solution in different regimes, 


it is useful to define (see also Eq. (5.19)): 

1 

\" 1 




^/,o 


U.N 


y(l - -fcfVmJ 
<5j^(|l-7c|5p)“' = 


7c 


U.N 


plume limit parameter (5.66) 

jet limit parameter (5.67) 


which are comparing Ufj, = Uo/(j) with Ug = g/N 925 m/s and 7 c with 1. As we 
will show in the next section, when is small (f/^ -C Ug and 7c < 1) the solution 
has mainly a plume-like behavior, on the contrary, when (5j -C 1, the solution 
behaves manly as a jet. 


®Here r(a:) is the Gamma function. 
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When we are in the plnme limit regime ((5p -C 1), any power of Ur, 
simplified to (see Eq. (5.59)): 

= |1 - (l + ^9’ = |1 - (l + iS'l + 0(^9) . 


can be 

(5.68) 


This approximation, leads to the limit 


<?max ^ 


^^m(l - 7c) 




^max - 5p 


-1 


f ~ 

J min — 


' 27 c - 1 if 7c < 1 
if 7c > 1 • 


if 7c < 1 
if 7c > 1 


(5.69a) 

(5.69b) 

(5.69c) 


Thus, in this regime we recognize two distinct behaviors: when 7c > 1 the multiphase 
plume is too heavy and slow to reach its height of positive buoyancy and it collapses. 
On the contrary, when 7 c < 1, the plume is able to reach its buoyancy reversal 
height and it can rise into the atmosphere. During its ascent, / varies approximately 
in [ 27 c — 1,1], while q and m reach a much larger value the more (5p is small. 

On the other hand, in the jet limit regime (5j -C 1) we have: 


ui = ((i-77 + 79’-W 

(5.70a) 

^ 1 + :^5i(i)iSj 

(5.70b) 

rUmax ~ 1 + - (1 - 

(5.70c) 

f ■ ~ 

jmin — 

(5.70d) 


In this case q and m reach maximum values near 1, while / decreases the more the 
more is small. 


5.6.1 Plume height 


Eq. (5.65) gives us the opportunity to write an analytic expression for the 
maximum height reached by a plume described by Eqs. ( |5.56[ ). Indeed, the maximum 
plume height (m=0) is reached when when / = /min (cf. Eq. (5.61)). Thus, by 
substituting / = /min in the integral lower limit, and performing a change of 
variable in the integral with x = {f — ^c)/V^m, we obtain (see dehnition for IXm in 
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Fig. 5.4 : Contour plot of the plume height function \){xQ^qQ) defined in Eq. ( |5.71[ |. This function assumes its 
maximum in = Fi ~ 2.572, and it is a strictly decreasing function of go- When xq 1~ we are 

in the plume regime] a::o ^ 0 


Eq. (5.59)): 



(5.71a) 

(5.71b) 

(5.71c) 

(5.71d) 


where l)(a;o,go) is a function defined in [—1,1] x [0, cxo). It is worth noting that 
with this substitution the neutral buoyancy level height can be easily obtained by 
substituting the lower bound of the integral x = —1 with x = 0 (cf. Eqs. (5.60) 
and ([5)6lt). 


In Fig. 5.4 we represent the values assumed by f)(a;o, %) in (xq, %) e (“IT) x 
(0,1). We notice that this function has a maximum in P)(1,0) = Ei ~ 2.572. 
Approaching this point, the function increases suddenly. This figure must be read 
keeping in mind four main regimes: 1) a:o 1“ when 7c < 1 and T 1- this 
case we are in the plume regime near the singular point (a:o) % = (h 0), thus the 
column initially has enough momentum to reach its buoyancy reversal height and 
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enough enthalpy to rise until its maximum; 2) when 7c > 1 and (5p -C 1, we are in 
the collapsing plume regime near the point (xq, go) = (—1, 0); 3) when -C 1 we are 
in the jet regime, near the line Xq = 0. In general, 7c is the parameter controlling 
the column stability; when 7c < 1 then 0 < Xq < 1, the column is not collapsing 
and when xq —t 1 the column behaves as a plume, while xq —)■ O’*", the column 
behaves as a jet. 

The expression for the plume height we have found is the multiphase version 
of to that found in Morton HD]. The behavior of f) near (xo,go) = (bO) is the 
more interesting from a volcanological point of view, and it can be studied by 
using asymptotic expansion techniques for f (plume regime). In this case, 

Eqs. (5.71) can be highly simplihed. Indeed by using Eq. (5.66), we have: 

1 


xo = sign(l - 7c) (^1 - -6l + 0((5j)^ 1 --6l 

% = ii - 7ci +0(5^/^) - (1 - 7c) 5:^5^ = 


2v, P 

Vr. 


2v, P 


ap = (l-7c)-, 


see 


footnot^H 


(5.72) 

(5.73) 

(5.74) 


because 7c < 1 near xq = 1 . Moreover, if x ~ 1, the hypergeometric function can 
be approximated as follows: 


p p o9/4 

x^{x‘^) = y (1 — x^)5 dx ~ 2^/"^ y (1 — x)3 =-^(1 — x)3 + 5(1). (5.75) 


With these information and Op small enough, say 


< 2^/^ + ~ 2.5 , 


(5.76) 


it is possible to show that: 


3;o ^ 

J dx go + a:o5q(xo)-a:5g(a;^) " - Ti 1 - Ea (l + 5^/“^ , (5.77) 


71' "" 


-1 


where 


El = ^ / dx [ 5 ,( 1 ) - x5g(x^) 


~ 2.572 


-1 


Eo ~ 0.3802. 


In this “plume regime”, the analytic formulation for the plume height given in (5.71) 
simplihes to the hrst order approximation: 


^(1) ^ ^(1) /£ ^ 
Smax max/ 


El 


1 - Ea 1 + a 


Vq ap 6p 




(5.78) 


^Recall that Op = W, see Eq. (5.401. 
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O 




0.01 


6p = 0.001°'^ = 0.03 

6p = 0.1 


6p = 0 .1 


0.5 

( 1 )/, 


0.3 


^max '^max 


( 0 ) 


a 


P 


Fig. 5.5: 


Comparison of the exact formula Eq. ( |5.71[ l for the plume height of model ( |5.56[ l with the first order 
approximation Eq. ( |5.78| | over the zeroth order approximation Eq. | |5.81| >. 


while the zeroth order approximation is: 


(‘(0) = = 
Smax max/ 


Ti 


Vq 




(5.79) 


This last approximation holds in the limit hp —)■ 0, which is equivalent to the pure 
plume solution with initial mass and momentum equal to zero and finite initial fiux 
of buoyancy. 

In Fig. 5.5 we show the good behavior of Eq. (5.78) when hp < 0.3 and Op < 5. It 


i worth noting from Tab. 5^ that this parameter range is the most interesting from 
the point of view of volcanic plumes. Fig. |5.5| compares the first order, the zeroth 
order and the exact solution (5.71). It shows that the first order approximation 


behaves very well in the selected parameter range. On the other hand, we point out 
that considering the first order approximation instead of the zeroth order allows to 
avoid an error np to 100% when hp ~ 0.3 and Op = 5 (i)fmax — -^mix — O-^-^mix)' 


5.4 


since go oc 


We obse rve also that Fig. |5.5| is a zoom on the singnlarity at the bottom right of 
Fig. 


In the literature, the problem of obtaining the maximum plume height starting 
from the monophase (yc = 0) formulation of the plume model in a stratified envi¬ 
ronment, Eq (5.57) has been studied in Morton et al. |TT]. He found Cmax,M — 2.805 


in his non-dimensionalization. We can recover the same result in the zero order ap¬ 
proximation, by noting that the conversion factor from our non-dimensionalization 
to that used by Morton et al. m is 


Cm = = 28Vqap6p 


(5.80) 


from which Cmax,M = 2® Fi ~ 2.805. Turning to dimensional variables, at the zeroth 
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order we have recovered the famous relationship: 

^(0) = (^9Qo\^ 

Nl ) v^UoiVoV 


(5.81) 


telling that the maximum plume height to the power four is proportional to the 
mass flow rate times the enthalpy anomaly and inversely proportional to the cube of 
the Brunt-Vaisalla frequency. In the monophase case, when the Ricou and Spalding 
[13] entrainment model can be considered a good approximation for the dynamics of 
the first part of the plume, this result is valid even if the Boussinesq approximation 
is not valid (see Eq. ( 5.52[ )). 

In volcanological applications the zero order formula is widely used. We have 
found a correction to that formula, for the multiphase case in both the zeroth 
and hrst order formulation. In dimensional variables, the multiphase hrst order 
formulation of the plume height reads: 


1 






1-W 


l+(^] 

l2x[/2 i 


An 

12 


UpNoy 

. (t>*9 J 


(/)* = (!- 7c)0 = 0 - [XsEs,o + (Xe - '4’e)Yefl] . 


(5.82) 

(5.83) 


which strongly increase the accuracy of the plume height, keeping a simple analytic 
formulation. The only difference between the monophase and the multiphase 
formulation is in the factor (1 — 7c), through the substitution 0 —)■ 0*. 

We remind that this Taylor series approximation holds when <C 1 which is 
equivalent to Uo/<f> < q/Nq ~ 925 m/s. This last condition give us a lower limit for 
0 and than to the vent temperature: 


> 


UoNo 


9 


ATo ^ 


T, 


Q ;,0 


9 


(5.84) 


If the vent temperature is much smaller than this lower bound, than the plume 


behaves more likely to a jet, and integral (5.71) must be evaluated without the 
approximation hp 1. 

When we are in the opposite condition hj = —)■ 0 (jet limit), we have 

xo —)■ 0 <C 1. In this regime, the function l)(a;o,^o) does not have a strong 
singularity as in the case xq —)■ 1 (cf. Fig. 5.4) and Eq. (5.71) can be safely 


approximated at the zeroth order as (use the fact that x^{x^) ~ x in x G [—1, 0]): 


H ~ 


'^q ( + 90 + % 


go = 


T/,0. 


2n„ 


4iVo 

4xf/n 


(5.85) 


If also go ^ 1 this expression further simplifies giving the following expression for 
the maximum jet height: 




xNq) 


(5.86) 


As a hrst order approximation one can use io — bo and invert this expression to 
hnd the inlet velocity from the jet height. 
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5.6.2 Neutral buoyancy level and plume height inversion 

By recalling that the neutral buoyancy level (nbl) is reached when / = 0, it is 
easy to modify Eqs. (5.71) and (5.78) to hnd ifnbr 


-^0 1 


f)nbi(a:o, go) = 


Vq{VmVffi)i 

dx 


' VmQ- - IcY + ^/,o ' 
. ^/,0 


f)nbi(a:o, go) 


^S/4 = 


72 

Bi 


qO+Xodq{xl) -X^q{x^) 


0 


Vq Up 6p 


Bnbl ~ r2 ( 1 + ttp^ ) 5 


Bnbl = 1 — 


72 Bi 


dx 


dq{l) - Xdqix"^) 


~ 0.7596. 


(5.87) 

(5.88) 

(5.89) 

(5.90) 


-1 


Thus we have found a hrst-order modihcation of the result of Turner 

hL^I ^ ^ B2(1 - B„bi) 


H, 


( 1 ) 

nbl 


nbl 


rf.b, 


1 + Clp' 




(5.91) 


At the zeroth order we hnd i7®x/-^nbi = 1/Bnbi — 1.316 in agreement with 
Bfmax/TBnbi = 1.3 obtained by Turner [TH] . 

This result is telling us that the ratio between the maximum plume height and 
its neutral buoyancy level is a constant B^i — 1.3 when is small enough, and it 
grows with 5^^. 

The neutral buoyancy level of a plume can be observed by measuring the height 
where the plume umbrella begins to spread up. If we know TTnbi, 77max, -^o — and 
the entrainment Vq = 2x, it is possible to invert Eqs. (5.78) and (5.89) in order 
to hnd 6p and Op or equivalently Uq, 0 and /do- Dehning hnbi = -f^maxZ-f^nbi and 
^max we hnd 


(cip) ^ + (Up) — cih 

(^nblBnbl - 1) 
BlB2hnbl(l — Bnbl) 


— 


(flp) " - ■: 


a 


h 


1 - 0.41a? + 1.4a? + l.OOa^ + a? 


Oh -r 

Bihnbl(l — Bnbl) 

5p = ^-7-7 (a 


Ua = 


nghniax(hnbl 1) 

ioNo 


VqQ,p6p 


7 = (1 - 7c)0 = 




(5.92a) 

(5.92b) 

(5.92c) 

(5.92d) 

(5.92e) 

(5.92f) 


VqQ ap^p 

a well posed problem when hnbi > B^i — 1.316. The hrst equation can be so lved 
looking for the unique positive root with respect x = (ap)“^/^ (cf. Eig. 5.6). In 
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Fig. 5.6: Root of Eq. (|5.92a^ as a function of aj, and its analytic approximation, Eq. l|5.92c[|. 


Eq. (5.92c) we give an approximate analytic solntion which has a good behavior 


both in the asymptotic {ah —)■ 0 and ah oo) and intermediate regime (0.5 < 
ah < 5). In conclusion, the first order approximation for the plume height gives an 
additional information allowing to hnd both Uq and (p* in contrast with the zero 
order approximation which needs an additional hypothesis on 0* to give the mass 
flux. 

In order to fix ideas, we give an example. Suppose to have a monophase air 
plume with Uq = 30 m/s, Tq = 373 K (00 = 0.947 kg/m^), bo = 0.1 m ejected 
in an atmosphere with Tq, o = 300 K, po = 101325 Pa and iVo = 1.015 * 10“^ 
Hz. Solving Eqs. (5.13) with the Ricou and Spalding [H] model (x = 0.14), we 
obtain — 1387.2 and = 1-3461, slightly bigger than Pjj/j ~ 1.316. 


Now, substituting i7max/^o, ^nbi, Vq = 0.28 and io = 0.1 m in Eqs. (5.92), we can 


invert the problem recovering the initial velocity and density. With our first order 
approximation, we obtain: 


Rq, inverted — 28 m/S 

00,inverted ^ 0.887 kg/m^ , 

with less than 10 % of error with respect to the “real” values. 


(5.93) 

(5.94) 


5.7 Analytic solution for a non-Boussinesq plume in a 
stratified environment 

In this section we want to find an analytic solution approximating the behavior 


of model (5.1) in its complete form, from the vent elevation up to the neutral 


buoyancy level. The strategy that we will follow here will bring to an update of the 
results we have presented in Cerminara et al. [3]. 


Both Eqs. (5.26) and (5.41) admit the same asymptotic solution fulhlling the 
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initial condition g(0) = 1 


?(C) = 




where Og = 


4Vn 


5vm(l - 7c) 


(5.95) 


Thus this solution approximate the plume model (5.1) in both the Boussinesq and 


non-Boussinesq regime. The difference between these two regimes appears in the 
asymptotic solution when we choose which first integral of motion to use, either U 


(Eq. ( |5.38D ) or Urs (Eq- ( |5.25[ )), thus in the form of m: 
1 


m(C) = 


Qq 


- 1 ) +1 


or 


^(0 = \ — [(^c(g(C)) - ^c(i)] +1 


with 


Uq) = q^- [q - qx in(|g + qx\)] 

/c 


(5.96) 

(5.97) 

(5.98) 


These asymptotic expansions are equivalent to Eqs. (5.46), with correct initial 


conditions m(0) = 1 and g(0) = 1. In what follows, we will use the latter Eq. (5.97) 
as asymptotic expansion for the momentum flux, because it works better than the 
former equation in the non-Boussinesq regime. Indeed, even if this solution has 


been found by applying the approximation g S> 1 to Eqs. (5.1), we want to extend 
its applicability to plumes in non-Boussinesq regime. We will describe a strategy to 
hold this task, after having introduced atmospheric stratification. 

The only difference between Eqs. (5.35) - from where we have extracted the 
latter asymptotic solution - and the Eqs. (5.56) - for a stratihed atmosphere - is the 


variability of /(C)- hr the former system / is considered as constant and equal to 1, 
while in the latter one it is considered as a function / = /(C)- However, we have seen 
in the previous section that f{z) is a slowly varying function, because Vf^ is usually 
very small with respect to the rate of variation of the other equations involved, 
namely Vg and Vm- Thus, one strategy to look for an analytic solution of the problem 


in a stratihed atmosphere could be to consider the asymptotic solution (5.95) valid 


also for problem (5.56), and use it for hnding /(C)- In particular, substituting q{Q 


in (5.56c), we obtain: 


/(C) = 1 - 


^/,o 


2(1 - Jc)Vr 




(5.99) 


in Eq. (5.59) 


with m(C) defined in Eqs. (5.96). Now, we recall the first integral of motion found 


Um = (1- Icf + - = (/- Icf + — 

Vm Vm 


(5.100) 


In Eqs. (5.461 are the asymptotic solution of system (5.35), written in a form such that it is 


possible to find the virtual radius b^. However, that solution does not fulfill initial conditions for 
q and TO. To write an asymptotic solution respecting the initial condition it is more convenient to 
use q{C) in the form given in this section. 
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parameter 

[forced Plume] 

[Santiaguito] 

[weakPlume] 

[strongPlume] 

Cmax 

1665 

23.98 

160.6 

20.68 

y(asy) 

Smax 

1487 

21.79 

139.8 

19.65 

Cnbl 

1264 

18.36 

118.6 

13.58 

^(asy) 

Snbl 

1145 

16.53 

106.1 

14.55 


Tab. 5.5: The main parameters defined in this section for the four plume examples of this thesis. 


and we try to substitute Eq. (5.99) in it. We find: 


Wo, 


Vffl2 




2\2 


(1 — m ) 


This result differs from Eq. (5.59) just because of the term 


(5.101) 


-'/.o 


4(1 - 7c) 


2y2 




2\2 


(5.102) 


where we have used the definition of hp = u/,o/(l — The latter term is 0(6^), 

thus it can be disregarded in the plume regime (5p -C 1) with respect the other two 
terms in the right-hand-side of Eq. (5.101), which are respectively 0(1) and 0((5p). 
By noting that Um is approximatively conserved by the asymptotic solution found 
in this section, we have corroborated the fact that this solution is approximating 
the complete solution in the plume regime. 

Having the enthalpy flux evolution /(C), it is possible to calculate the maximum 
plume height and neutral buoyancy level by using mmax and /min given respectively in 


Eqs. (5.60) and (5.61). In Tab. 5.5 we recall the maximum plume height and neutral 


buoyancy level as obtained from model (5.1), comparing it with the asymptotic 

-rpciiifc 

lesuns . 

Now we move to face the non-Boussinesq regime. The strategy we proposed 
in Cerminara et ah [5] is to use the asymptotic solution in the complete inversion 
formulas for U, b, (3, T^, W and W reported in Eq. (5.5). The behavior of this 
approximation is showed in Eigs. 5/7, 5.8 IffOllSTOl There we notice that the 


solution works surprisingly well for all the presented plumes. In particular, the 
temperature and density profiles are well captured for all the cases. The best 
behavior is recorded in the non-Boussinesq monophase plume (recall 0 = 0.893). 
The asymptotic solution behaves worse for the plume radius and the plume axial 
velocity in the upper part, where the stratification play the most important role. 
Anyway, the plume height is captured with less than 10 % of error for all the plumes. 


Systematically, the asymptotic mass flux is overestimated with respect model (5.1). 


This error present with more evidence in strongPlume, and directly reflects in the 
underestimation of the mass fractions along the plume axis. 
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I - <1 - m - / - q - m . f\ 

(a) mass, momentum and enthalpy fluxes 






Fig. 5.7: [forcedPlume]: Vertical evolution of the non-dimensional fluxes m, f (log-log scale), of the plume 
radius b (log-log scale) and of the dimensional physical parameters U, /?, T^, Ve(s)? (linear-log) scale. 


Solid lines correspond to the numerical solution of model 
the analytic asymptotic solution Eqs. ( |5.95[ l, ( |5.97[ l, |5 


el 

■99[|. 


while dashed lines are evaluated by using 
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I - <1 - m - / - q - m . f\ 

(a) mass, momentum and enthalpy fluxes 






Fig. 5.8: [Santiaguito]: Vertical evolution of the non-dimensional fluxes m, f (log-linear scale) and of the 
dimensional physical parameters U, 6, /3, T^, Solid lines correspond to the numerical solu¬ 
tion of model while dashed lines are evaluated by using the analytic asymptotic solution 

Eqs. ( |5.97[ l, |5.99| |. 
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I - <1 - m - / - q - m . /I 

(a) mass, momentum and enthalpy fluxes 





P [kg/m3] 


(d) plume density 


z [km] 



Fig. 5.9: [weakPlume]: Vertical evolution of the non-dimensional fluxes g, m, / (log-linear scale) and of the 
dimensional physical parameters U, 6, /3, T^, Solid lines correspond to the numerical solu¬ 
tion of model while dashed lines are evaluated by using the analytic asymptotic solution 

Eqs. ( |5.97[ l, |5.99| |. 
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(a) mass, momentum and enthalpy fluxes 





(d) plume density 




(e) plume temperature (f) plume mass fractions 


Fig. 5.10: [strongPlume]: Vertical evolution of the non-dimensional fluxes g-, m, / (log-linear scale) and of 
the dimensional physical parameters C/, b, /3, T^, Ve(s)' Solid lines correspond to the numerical so¬ 
lution of model jgp while dashed lines are evaluated by using the analytic asymptotic solution 
Eqs. l |5l97| l, ( |5.99| |. 
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6 Comparison between results of 3D and integral plume 
models 


Integral models for plumes describe the evolution with height (the axial unity 
vector being z) of three main variables: the flux of mass, momentum and buoyancy. 
The purpose of these kind of models is to reproduce - as accurately as possible 
- the behavior of these three parameters under the hypothesis that the plume is 
stationary. Moving to the 3D models, they give us the plume variables as a function 
of time and space. In order to compare results, we have hrst of all to average the 
3D result over a time window where the solution can be considered stationary. The 
second step to do in order to coherently compare the two kind of models is to define 
the three fluxes also in the 3D case. We choose to dehne it as described below. 

Given O x T, the space-time domain, we hrst average over T a generic 3D 
variable 'ijj{x,t): 

( 6 . 1 ) 

For keeping the notation as simple as possible, in this section we use (■) in place of 
(•)j. We dehne a plume subset f2piin(z) C Dz, where is the plane orthogonal to 
i at height 2 ;. Subset Dpim is identihed by two thresholds: the averaged mixture 
velocity has positive axial component and the mass fraction of a tracer ^tracer is 
larger than a minimum threshold Umi-a'- 


f^plm — {{Xl) G D 2 ; I Itm ' Z ^0 and I/tracer ^ 2/min} • (6-2) 

We refer to the integral over this domain as: 


V’(^) = ('0(a;))npi,n = / dxidx2'il^{xi,X2,z). (6.3) 

In particular we dehne respectively the mass hux, the kth mass fraction, the 
momentum hux and the buoyancy hux as follows: 

ttQ = {pm = 7r/?f/6^ (6.4a) 

T^Qk = {PmVk = 'xldYkUb'^ (6.4b) 

TtM = (pm{Um ■ (6.4c) 

' ' ^^plm 


'kF 


1 1 + J2k{Xk 
\1 + Hkii^k 


^)yk 

^)yk 




= TT 


^plm 


l + K 


-a 


P 


/3 j Ub^ , 

(6.4d) 


where = J^kii^k — 1)W , Y^ = J^kiXk — 1)W and fc G J U 3 (with nil gas constant 
of the solid phase -ipj = 0). Moreover, a{z) = {pa{x))Q^^^. We choose this method 
for obtaining the one-dimensional integral huxes because of two reasons: 1) it is 
the three-dimensional counterpart of what we have dehned in Secs. |^and|^ thus 
it holds even in non-Boussinesq regime 2) it is independent on the shape of the 
radial prohle of the plume. 


similar approach for the Boussinesq regime has been developed in Kaminski et al. [7]. 
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By defining = Y^Q and = Y^Q, we can recover the plume variables by 
using the same inversion formulas given in |5.5[ We recall them in their dimensional 
form: 


• plume radius b{z) = 



• plume density /3{z) = a 

• kth averaged mass fractions Yk{z) = % 


• plume temperature T{z) = Ta q^q 


• plume velocity U{z) = ^ 

• entrainment coefficient >c{z) = 

where (■)' is the derivative along the plume axis and Tq = p/Rad is the atmospheric 
temperature profile. 

It is worth noting that the methodology described in this section allows plume 
modelers to coherently compare results obtained from one-dimensional integral mod¬ 
els with data obtained from complex three-dimensional simulations. Moreover, the 
entrainment coefficient x - the key empirical parameter for one-dimensional models 
- can be easily obtained for three-dimensional helds. In |2] we give some example of 
the results we obtain by using this averaging procedure for the post-processing of 
three-dimensional plume simulations. We have used the same procedure also for 
the lAVCEI (International Association of Volcanology and Geochemistry of the 
Earth Interior) plume model intercomparison initiative |1], consisting in performing 
a set of simulations using a standard set of input parameters so that independent 
results could be meaningfully compared and evaluated, discuss different approaches, 
and identify crucial issues of state of the art of models. 
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A Notation 


a acceleration 
b plume radius 
c speed of sound 
C specific heat 
Cd drag coefficient 
Cp specific heat at constant pressure 
Cv specific heat at constant volume 
C compressibility of the velocity 
field: (|V-t.nj,/(|Vtt|2)^, 
d particle diameter 
d spatial dimension 
D vent diameter 
T) strain rate tensor 
e internal energy per unity of mass 
B total energy per unity of mass 
£ kinetic energy per unity of mass 
spectrum 

fj drag force per unity of volume 
acting on the jth particle class 
F buoyancy flux 

2 -Fi, Gauss hypergeometric functions 
g gravitational acceleration norm 
g' reduced gravity 
g gravitational acceleration vector 
g gravitational acceleration versor 
“K enstrophy per unity of mass 
h enthalpy per unity of mass 
-f^max volcanic plume maximum height 
hfnbi volcanic plume neutral buoyancy 
level 

i index running over all the chemi¬ 
cal components in the fluid phase 
I number of chemical components 
in the fluid phase 
J set of all the indexes i 
I identity tensor 

j index running over all the parti¬ 
cle classes 

J number of particle classes 
3 set of all the indexes j 
k wavenumber 


kg thermal conductivity 
K kinetic energy per unity of mass 
Kt subgrid-scale kinetic energy per 
unity of mass 
L length scale 
m mass 

N number of grid cells 
N Brunt-Vaisalla frequency 
p pressure of the fluid phase 
q heat flux 
r radial coordinate 
f radial unity vector 
R gas constant 
Q mass flow rate 

Qj heat per unity of volume ex¬ 
changed from the fluid phase to 
the jth particle class 
Qw release of thermal energy from 
the vent 

Q subgrid-scale diffusivity vector 
for the temperature 
S source term 
S rate-of-shear tensor 
S vorticity tensor 
t time 

T temperature 
T stress tensor 
7 temporal domain 
u velocity vector 
U velocity scale or mean plume ve¬ 
locity 

f/e entrainment velocity 
V volume 

w particle settling terminal velocity 
W WALE subgrid model operator 
X position vector 
y mass fraction 

y subgrid-scale diffusivity vector 
for the mass fraction 
axial coordinate 
i axial unity vector 

a density of the atmosphere 
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/9 gas-particle mixture density for 
the integral plume model 
/3p density ratio parameter 
7 adiabatic index of the gas mix¬ 
ture 

7 c stability of the plume column 
6 grid scale 


Nu Nusselt number 
Pr Prandtl number 
Prt subgrid-scale turbulent Prandtl 
number 

Re Reynolds number 
Ri Richardson number 
St Stokes number 


Ax smallest space scale of the dy¬ 
namical problem 
e volumetric concentration 
et subgrid-scale energy dissipation 
( non-dimensional axial coordinate 
7 k Kolmogorov length scale 
7 x entrainment function 
6 atmospheric thermal gradient 
azimuth angle 

K dispersed on carrier mass ratio 
X entrainment coefficient 
At Taylor microscale 
u fluid kinematic viscosity 
^ smallest resolved LES length 
scale 

fx fluid dynamic viscosity 
/ib fluid bulk viscosity 
/it subgrid-scale eddy viscosity 
p bulk density 
p density 
g density scale 
r typical time scale 
Te eddy turnover time 
Tp Kolmogorov time scale 


(■) cell faces averaging 
{■)q space domain averaging 
(■)t temporal domain averaging 
{■)j jth mass fraction weight average 
over the domain 
(■) filtered quantity 
(■) Favre-filtered quantity 
(■)dg dusty gas 
(■)e ejected gas phase 
(■)f fluid phase 
(■)g gas phase 

(■)i ith chemical component of the 
fluid mixture 
(■)j jth particle class 
(■)r correction due to particle decou¬ 
pling 

(■)rms root mean square 
(■)s solid phase 
(■)sth Sutherland law 
(■)q, atmospheric 
(■)m gas - particle mixture 
(•)/3 gas - particle mixture (integral 
model) 


V molar fraction 


(pc drag correction function 
X ratio between specific heats 
tjj ratio between the gas constants; 

generic function 
n spatial domain 


Co Courant number 
Ec Eckert number 
Eu Euler number 
Fr Fronde number 
Ma Mach number 
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B Gauss hypergeometric functions 


Gauss hypergeometric functions 2 -^i([‘, ■] i are useful in order to perform 
integrals of the form; 


J (x^ — a)^ dx . 

2 -Fi([-, ■]; [■]; x) is the hypergeometric fnnction dehned when a; < 1 as: 

{a)n{h)nX^ 


2 Fi(a, 6;c;a;) = ^ 

.1 

[G>]n. 


n=0 vGn rt. 


n = 0 


a{a + 1 )... (a + n + 1) n > 0 . 

In thesis we have to deal with integrals in which c = 2, thus we dehne 

db{x) = 2-^1 


- 6 ,- 

2 . 


■,x 


(5b{x) = 2 F 1 

so that 


A- 6 -- 




; X 


X 


a — x^ydx = a°x^b — ] + C 
\ a / 


x‘^ — a^dx = 


X 


l+2b 




1 + 26 \x 


+ C 


if x"^ < a 
if x‘^ > a . 


(B.l) 


(B.2) 

(B.3) 


(B.4) 

(B.5) 

(B.6) 

(B.7) 


It is worth noting that and C5fe(l) are hnite and them value is tied to the 

Gamma function r(a;) as: 


yirr(l-6) 

' 2r(3/2-6) 

(B.8) 

, 22''0Fr(l-26) 

’ r(l/ 2 - 26 ) ■ 

(B.9) 
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